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^ . Abstract 

I— 1 1 Waves scattered by a weakly inhomogeneous random medium contain a predominant single 
rlln . 

^p^' scattering contribution as well as a multiple scattering contribution which is usually neglected, 
c/2 ! 

^ ' especially for imaging purposes. A method based on random matrix theory is proposed to sep- 

CJ , arate the single and multiple scattering contributions. The experimental set up uses an array of 
c/j ■ 

O \ sources/receivers placed in front of the medium. The impulse responses between every couple of 
'c/3 ■ 

transducers are measured and form a matrix. Single-scattering contributions are shown to ex- 

^ : 

Q-i. hibit a deterministic coherence along the antidiagonals of the array response matrix, whatever the 

\ distribution of inhomogeneities. This property is taken advantage of to discriminate single from 

^T) ' multiple-scattered waves. This allows one to evaluate the absorption losses and the scattering losses 

0^ ■ separately, by comparing the multiple scattering intensity with a radiative transfer model. More- 

, over, the relative contribution of multiple scattering in the backscattered wave can be estimated, 

o ■ 

I which serves as a validity test for the Born approximation. Experimental results are presented 

T— I ! 

^ ' with ultrasonic waves in the MHz range, on a synthetic sample (agar-gelatine gel) as well as on 

^ . breast tissues. Interestingly, the multiple scattering contribution is found to be far from negligible 

H ; 

5r ! i'^ the breast around 4.3 MHz. 
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I. INTRODUCTION 

Standard imaging techniques, such as ultrasonic echography [li, radar or optical co- 
herence tomography [sl, are based on the same principle. One or several source (s) emit(s) 
a wave into the medium to be imaged. It is reflected by the inhomogeneities of the medium 
and the backscattered wave is measured by the same or other sensor (s). It contains two 
contributions: 

• A single scattering contribution (path s in Fig{T]): the incident wave undergoes only one 
scattering event before coming back to the sensor (s). This is the contribution which is used 
in imaging because there is a direct relation between the arrival time t of the echo and the 
distance d between the sensors and the scatterer, t = 2d/c (c is the sound velocity). An 
image of the medium reflectivity can be built from measured signals. 

• A multiple scattering contribution (path m in Fig{T]): the wave undergoes several scattering 
events before reaching the sensor. Multiple scattering occurs when scatterers are strongly 
diffusive and/or highly concentrated. There is no correspondence between the arrival time 
t and the position of a scatterer. Thus, classical imaging fails in multiple scattering media 

Standard imaging techniques rely on the single-scattering assumption (flrst Born approx- 
imation). However, there is no such thing as a purely single scattering medium. A multiple 
scattering contribution always exists, albeit negligible compared to single scattering. Nat- 
urally for imaging purposes, one tries to reduce the influence of multiple scattering, for 
instance by choosing an appropriate frequency domain where multiple scattering is not too 
strong 8|. Focused beamforming with an array of transducers, or more generally synthetic 
aperture techniques are also a way to enhance the single scattering contribution. It 
should be noted that even though multiple scattering is considered as the enemy of classi- 
cal imaging techniques, studying it may bring additional information about the scattering 
structure. Indeed, a wave undergoing multiple scattering can be thought of as a random 
walker j^, with two essential parameters: the elastic mean-free path and the diffusion 
constant D. Measuring the se p arameters is a way to characterize the microarchitecture of 
the scattering medium 10|, . Yet in weakly inhomogeneous media where the flrst Born 
approximation is reasonably valid (especially human soft tissues probed by ultrasound in 
the MHz range), it is a challenge to study multiple scattering parameters because of the 
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predominance of single scattering. 

Recently, an original technique has been proposed to separate the single-scattered echo 
of a target drowned in a predominant multiple scattering background jsl.]]^. The method 
was based on a matrix approach. It has been successfully applied to target detection and 
imaging in highly scattering media. In this paper, we are also interested in discriminating 
single-scattering and multiple-scattering contributions from the total response of an un- 
known medium, based on matrix properties. However, the present approach is different 
from earlier works 0, Q, both in terms of method and applications. The situation we 
consider here is exactly the opposite: we want to extract the multiple scattering contribu- 
tion from predominantly single-scattered waves in a weakly scattering media. Moreover, the 
method is based on a singular value decomposition applied to the antidiagonals of the array 
response matrix, and not to the array response matrix itself. The distinction between sing. 



and multiple scattering subspaces is then performed using random matrix theory 13|, Il4j . 
as it will be detailed in the next sections. 

The interest of this work is twofold. First, once single and multiple scattering contri- 
butions are isolated, the proportion of multiple scattering within the wave response of the 
medium can be evaluated. This figure can be used as an indicator for the validity of the 
single-scattering (Born) approximation, which is the basis of classical imaging techniques. 
Note however that the purpose of this study is not to improve imaging of weakly scattering 
media. The second interest of this work is to provide a new tool for the characterization 
of weakly scattering media. More precisely, we will show how the multiple scattering con- 
tribution, once it is isolated, can be taken advantage of in order to estimate the scattering 
mean-free path 1^ independently from the absorption mean-free path /„, thus discriminating 
absorption and scattering losses. 

The experimental results we present were obtained with pulsed ultrasonic waves firstly 
in a synthetic medium (agar-gelatine gel) around 3 MHz then in breast tissues around 4.3 
MHz, but the principle of the technique can be applied to all fields of wave physics {e.g 
seismology, electromagnetism, acoustics etc.) for which the multi-element array technology 
is available and provides time-resolved measurements of the amplitude and the phase of the 
wave field. 
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FIG. 1: Experimental setup: a 125-element linear array is placed in front of a random medium at 
a distance a. The whole setup is immersed in a water tank. 

II. TRANSDUCERS' ARRAY CONFIGURATION 



We use a N-element ultrasonic array (here, N = 125). The array is placed at a distance 
a from the random scattering sample under investigation (see Fig. [I]). The first step of the 
experiment consists in measuring the inter-element matrix. A sinusoidal burst of length 
St, at the central frequency fc, is emitted from transducer i into the scattering medium. 
Typical values here are 5t ~ 1 /is and fc ^ "i MHz. The backscattered wave is recorded 
with the transducers of the same array, which yields a set of impulse responses hij{t) 
{j = 1,...,N denotes the receiver index). The operation is repeated for the N emitting 
transducers. The responses hij{t) form the N x N impulse response matrix H(t). Because 
of reciprocity, hij{t) = hji{t) and H(t) is symmetric. A short-time Fourier analysis of the 
impulse response matrix H is performed. The time signals hij{t) are truncated into At-long 
overlapping windows : fey(T,t) = hij{T - t)WR{t) with Wnit) = 1 for t G [-At/2 , At/2], 
WFi{t) = anywhere else. The value of At is chosen so that signals associated with the same 
scattering event (s) within the medium arrive in the same time window 15|]. Typical values 
here are At ~ 10 fis. A Fourier analysis of K(T, t) is achieved by means of a discrete Fourier 
transform. A response matrix K(T, /) is finally obtained at each time T and frequency /. 
The single and multiple scattering contributions can now be discriminated with the help of 
a matrix manipulation. 
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III. THE SIGNATURE OF SINGLE SCATTERING 



When studying the array response matrices K(T, /), the predominance of single scat- 
tering manifests itself by the presence of a long-range deterministic coherence along the 
antidiagonals of the matrix, whatever the distribution of scatterers s], Q, 15|. As an exam- 
ple, Fig|2] shows the real part of one of the matrices K, in the case of a synthetic medium 
(Agar-gelatine gel) which is enough weakly scattering for the Born approximation to be valid. 
Even if the inhomogeneities are randomly distributed, K obviously exhibits some kind of 
coherence along its antidiagonals (i.e., for matrix elements kij such that i + j = constant). 
This coherence is a typical signature of single scattering, and it vanishes when multiple 
scattering dominates. This has been thoroughly explained in 
main argument here. 



15| , we briefly recall the 




J 

FIG. 2: Real part of matrix K obtained in a gel (5% gelatine, 3% agar-agar) at time T = 114 //s 
and frequency / =3.05 MHz. The source-sample distance was a = 50 mm. 



Generally, the kij{T, f) can be written as the sum of single and multiple scattering con- 
tributions: 

hATJ) = kf^{TJ) + kff{TJ) (1) 

Under the paraxial approximation, the distance between the origin (0, 0) and an observer 
(x, R) located slightly off-axis (x << R) is y/R"^ + x'^ = R + x^ /{2R). As a result, the phase 
shift undergone by a wave travelling from a source with coordinates (0,Xj), scattered by a 
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point {Xd,R) and received in the plane of the source at (0,Xj) reads : 



exp [2jkR] exp 



jk 



2R 



exp 



jk 



2R 



with k the wavenumber. The quadratic phase terms can be factorized since 



(xi - Xd) + (xj 



Xd 



[Xi 



+ 



[ OC 2 I 3C 



2Xd 



(2) 



Consider an ensemble of scatterers randomly distributed. As long as only the first and 
the last scattering of every scattering path are identical (which is naturally the case, if only 
single scattering takes place) the coefficients of the array response matrix at time T and 
frequency / will be proportional to: 



N, 
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(3) 



with R = cT/2 and A^^^ the number of scatterers contained in the isochronous volume. Ad 
depends on the reflectivity of the scatterer. Ad and Xd are random variables so kfj is itself 
random. Interestingly, applying the factorization of Eqj2] to Eqj3l a deterministic relation 
arises along the antidiagonals of K^: 
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(4) 



with 5x the array pitch and 2m6x denotes the distance between two array elements {i—m and 
i + m) on the same antidiagonal. EqJH implies that as long as there is only single scattering, 
there must be a form of coherence, a long-range deterministic relation, between the elements 
of the array response matrix, whatever the realisation of disorder. On the contrary, when 
multiple scattering occurs (except for recurrent scattering paths [16|, but this contribution 



is negligible in weakly scattering media), the elements k^J cannot be factorized, and there 



is no such long-range deterministic coherence 
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IV. SEPARATION OF SINGLE AND MULTIPLE SCATTERING 

The key to separate single (K^) and multiple (K'^) scattering contributions is the par- 
ticular coherence of along its antidiagonals. In previous works js], , was extracted 
from K by projecting the antidiagonals of K along the vector [f3m] of EqJH But the simple 
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form taken by EqJH results from a series of assumptions (paraxial approximation, pointlike 
array elements and scatterers) which do not all apply to our experimental configuration. In 
order to separate and K^, the method proposed in this paper is much less restrictive. 
We do not assume that EqS] exactly applies; we only assume that because of single scatter- 
ing there must be a deterministic coherence between the antidiagonal elements of K^, but 
we do not suppose we know its exact form. 

Under these conditions, the separation between and will essentially rely on a 
singular value decomposition (SVD) of the antidiagonals of K. This separation is a three- 
step process: 

• Rotation of each matrix K and construction of two sub-matrices Ai and A2. 

• Filtering of matrices Ar (r = 1,2). Ar is decomposed as the sum of two matrices: 
Ar = A^ + A^, where A^ and contain respectively the single and multiple scattering 
signals. 

• Construction from A^ and A^ of the single and multiple scattering matrices and K'^. 

The first and third steps (rotation of data) have already been presented in previous works 
[sl, [1^ and will be briefly recalled in Sec ilV Al and Sec ilV Cl On the contrary, the second 



step (SVD of antidiagona. 
the previous approach 
in Sec ilVBl 



^constitutes the core of the method and differs completely from 
12l | . The corresponding matrix operations are explained in details 



A. First step 

A rotation of matrix data is achieved as depicted in FigjSl It consists in building two 
matrices Ai and A2 from matrix K = [kij]: 

Al = [aiuv] of dimension (2M - 1) x (2M - 1), 

such that ai[u, v] = k[u + V — 1, V — u + 2M — 1] (5) 
A2 = [a2uv] of dimension (2M - 2) x (2M - 2), 

such that 02 [m, v] = k[u + v,v — u + 2M — 1] (6) 

where M = [N + 3) /A. Here = 125 and so M = 32 is an even number. The columns of 
matrices Ai and A2 correspond to the antidiagonals of K (see Figl3]). In the next subsection, 
we will no longer make the difference between matrices Ai and A2 because they are filtered 
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FIG. 3: Example of a matrix K of dimension N = 17. The black points represent the elements kij of 
K. The antidiagonals of K are the columns of matrices Ai and A2. Circles and squares represent 
the elements of Ai and A2, respectively. Once single and multiple scattering contributions are 
separated, the final matrices and have (2M — 1) x (2M — 1) elements (central square). 

in the same way. They will be called indifferently A. L is the dimension of A. For matrix 
Ai, L = 2M — 1; for matrix A2, L = 2M — 2. Because of spatial reciprocity, K is symmetric 
{kij = kji). Thus, A exhibits also a symmetry: each line of its upper part is identical to a 
line of its lower part. The symmetry axis is shown as a black line in Figj3]and corresponds to 
the diagonal of the matrix K. So, each column of the matrix A contains only M independent 
coefficients, even if its dimension L is larger than M. 



B. Second step 

A can be written as a sum of two matrices A^ and A^, which correspond respectively 
to the single and multiple scattering contributions 

A = A^ + A^ (7) 

Contrary to previous works jsl, Q, the technique we propose in this paper consists in sep- 
arating single and multiple scattering by achieving the singular value decomposition (SVD) 
of the matrix A. The SVD decomposes a matrix into two subspaces: a signal subspace (a 
matrix characterized by an important correlation between its lines and/or columns) and a 
noise subspace (a random matrix without any correlations between its entries). When the 
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SVD is applied to the matrix A, the signal subspace {i.e., the largest singular values) cor- 
responds to (the single scattering contribution characterized by a long-range correlation 
along its columns) and the noise subspace {i.e., the smallest singular values) corresponds to 
A^ (the multiple scattering contribution). 
The SVD of matrix A is given by 

L 

A = UAVt = ^AfeUfcVi (8) 

k=i 

U and V are square unitary matrices of dimension L. Their respective columns and 
correspond to the singular vectors associated to the singular value A^. A is a square diagonal 
matrix of dimension L, containing the real positive singular values in a decreasing order 
(Ai > A2 > ... > Xl)- Actually, A has only M non-zero singular values since it contains only 
M independent lines, hence EqjH] becomes: 

M 

A = UAVt = J]AfeU,Vl (9) 

k=l 

The issue is to determine which rank of singular value separates the signal subspace 
(single scattering) from the noise subspace (multiple scattering). If EqJUwere strictly true, 
the single scattering contribution A^ would be of rank 1 and only the first singular space 
associated to the first singular value Ai would correspond to the signal subspace. But when 
assumptions leading to EqJH do not strictly hold, A^ is no longer of rank 1 and several 
singular spaces (associated to the largest singular values) are needed to fully describe the 
signal subspace. This happens for instance when scatterers are not pointlike or when the 
paraxial approximation does not hold. We have to define a threshold to discriminate the 
signal and noise subspaces, with the help of random matrix theory (RMT) isl |l^ . 

By convention and for the sake of simplicity, the singular values A^ are first normalized 
by their quadratic mean 

A, = , (10) 

For a random matrix of dimension P x Q (with 1 << P < Q), whose entries are complex 
random variables, independently and identically distributed, the probability density function 
p(A) of the normalized singular values A^ is given by jl^ 



p(A) = ^J(A?,ax-A^)(A^-A^i,) (11) 



for Ajj^j]^ < A < Amax and otherwise, with 

-^max,min = 1 i a/P/ Q. (12) 

For random matrices of large dimensions, the singular value spectrum has a bounded support. 
In our case, A is a square matrix of dimension LxL. Yet, as it contains only M independent 
lines, it is equivalent to a rectangular M x L matrix. If it were truly random (which is 
expected to be the case of the multiple scattering contribution) its largest singular value 
should not exceed Amax = 1-71. This value is obtained from EqUll with P = M = 32 and 
g = L = 63. 

It should be noted that rigorously L and M are not large enough for the asymptotic law 
(Eq JTT]) to apply. Actually the first singular value Ai obeys a complicated law, known as 
the Tracy- Widom distribution jl3], which is of unbounded support. The probability for Ai 
to be larger than Amax can be computed: it is found to be ~0.08 for P = M = 32 and 
Q = L = 63. The presence of correlations between matrix entries also induce a deviation 
from EqJTT] as we will see later. For the sake of simplicity, we admit for now that within an 
acceptable probability of error, the singular values are upper bounded by Amax- 

According to EqlTJ A is the sum of a matrix A^ of rank p < M (associated to single 
scattering) and a matrix A'^ of rank M (associated to multiple scattering). Sengupta and 



Mitra [IJ] have shown that the (M—p) smallest singular values (linked to the noise subspace) 
exhibit the same distribution as singular values of a random matrix whose size is (M—p) x L. 
Let Amax denote the upper bound of the singular values distribution in the case of a random 
matrix of dimension (M — q) x L; we have: 

Ai^^ax = 1 + ViM-q)/L (13) 

From this property, one can propose a way to separate the signal and noise subspaces of 
A. We first consider the first singular value Ai upon normalization (Eq JTOj) . If Ai is larger 
than ASSax (EqMD, it means that the first singular space AiUiV[ is associated with the 
signal subspace. Then, we iterate the process and consider the second singular value. The 
Afc are once again renormalized, considering only the singular values from k = 2: 

A2 = , = (14) 
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The threshold value Amax to consider this time is the one obtained for a random matrix of 
size (M — 1) X L, i.e A^^^x (Eq JT3|) . If A2 > A^ax? the second singular space A2U2V2 is 
also linked to the single scattering contribution, and we iterate once more the process until 
rank p + 1 for which Ap+i < A^^x- Finally, we obtain a threshold rank p which allows to 
separate the signal (S) and noise (N) subspaces, 

p M 

S = ^ AfeUkVi ,1^= Yl ^fcUkVi (15) 

fc=i fc=p+i 

Ideally, S should be devoid of multiple scattering. This is not strictly true, because multiple 
scattering signals are not strictly orthogonal to the single scattering subspace. Let (t| and 
be the power of single and multiple scattering signals. In Appendix A, we show that 
the typical amplitude of the remaining multiple scattering contribution in S is au^/p/'^M 
(<< as). If we neglect this residual term, we have separated single and multiple scattering 
contributions: ~ S and A'^ ~ N. The whole separation process is summarized in FigJH 



SVDofA A = UylV 

9 = 1 



Determination 
of the superior 
bound of p(X) 



Normalization of 
singular values j 



1 



S4 



Comparison 



if 



Rank of s eparation : p-q-\ 
Single scattering: A^f» S = ^ \ 



Multiple scattering: A"f»N= 'Y^ A^U^V^ 



FIG. 4: Principle of the separation between the single and multiple scattering contributions. 

Note that a multiple scattering rate 7 can be directly measured from the singular values 
Afc of A. The sum of the square of all the singular values corresponds to the total intensity 
backscattered by the medium towards the transducers' array . Hence, a multiple scattering 
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rate 7 can be estimated from the singular values of A: 

7 = .0 (16) 



Until now, for simplicity we have implicitly assumed that along the antidiagonals of K'^ 
(or the columns of A'^) the matrix elements are completely decorrelated. However, experi- 
mentally short-range correlations may exist between elements, mostly because of mechanical 
coupling between neighboring transducers and of the coherence length of the diffuse wave- 
ield I12I. Il5|. Correlations between matrix elements can be taken into account theoretically 
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Consequently, the actual probability density function p(A) is more complicated than 

M 



the simple result of EqUH which modifies the upper bound A^ax [121, ll5|. In practice, A 



has to be computed numerically, based on a acceptable probability of error (typically 1%). 
The details are given in Appendix B. 

This technique of separation is based on the fact that the first singular value exceeds the 
value Amax, otherwise there is no separation between single and multiple scattering and the 
whole signal is considered to be associated with multiple scattering. So, this approach is not 
well suited for strongly diffusive media, i. e random media for which the multiple scattering 
contribution is predominant 

C. Third step 

The third step is the reverse of the first one. From A^ and A'^, two matrices and 
K^, of dimension (2M — 1) x (2M — 1), are built(see Fig|3]) with a change of coordinates, 
back to the original system: 

• if (z — j)/2 is an integer, 

then, = af'*'^[(^-j)/2 + M, {^+J)/2] 

• a {i — j)/2 is not an integer, 

then, k^'^'[t,j] = al'''[it - j - l)/2 + M, (^ + j - l)/2] 

contains the single scattering contribution (plus a residual multiple scattering contribu- 
tion) and contains the multiple scattering contribution. 
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V. CHARACTERIZATION OF A WEAKLY SCATTERING MEDIUM 



The experimental set up has aheady been described in SecJTTland is shown in Fig{Tl The 
experiment takes place in a water tank. The ultrasonic array has = 125 elements. The 
emitted signal is a sinusoidal burst of length 5t = 2.5 fis at the central frequency (3 MHz). 
The sampling frequency is 20 MHz. Each array element is 0.39 mm in size and the array 
pitch 6x is 0.417mm. The source-sample distance is a = 50 mm. The first random medium 
of interest is a gel composed of 5% of gelatine and 3% of agar. In this kind of medium and 



frequency range, the single scattering contribution is by far predominant [18| . The thickness 
L of the scattering slab is 100 mm. Once the inter-element matrix H is measured, the short- 
time Fourier analysis described in SecUl] yields a set of response matrices K(T, /). Then, 
the separation of single and multiple scattering is achieved as described in Sec JIVI 

Figj5] shows a typical experimental result, taken at time T = 114 and frequency / = 
3.05 MHz. Note that the separation rank between the signal and noise subspaces is here p = 
2, which confirms that EqJUdoesnot strictly hold.. exhibits the deterministic coherence 
along the antidiagonals (FigJSI^a)) which is characteristic of single scattering. Obviously, 
is very close to the raw matrix K (Fig|2]), since single scattering is predominant. As to K'^, 
it displays a random feature as expected for the multiple scattering contribution (Figj5]^b)). 
However, one cannot conclude that it originates in multiple scattered waves: it could also 
correspond to experimental noise. 




FIG. 5: Separation of single and multiple scattering contributions at time T = 114 /us and frequency 
/ = 3.05 MHz. (a) Real part of K^. (b) Real part of K^. 



In order to establish the multiple-scattering origin of K , we calculated the mean 



backscattered intensity as a function of the source-receiver distance X = Xj — Xi 



m5x 
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and the arrival time T 



(17) 



The symbol (.) denotes an average over the quantities in the subscript, here frequency and 
all source/receiver couples separated by the same distance X. In FigJHl the spatial 
dependence of /^^ is compared with the total intensity, at a given time T. Whereas the 
total intensity / shows no preferred direction, I^\X) exhibits a typical signature of multiple 
scattering: the coherent backscattering peak clearly arises around X = 0. This phenomenon 




-2 2 

X[mm] 



FIG. 6: The multiple scattering intensity I'^'^ (black dots) and the total intensity I (grey circles) 
are plotted versus X, at time T = 137 //s. The intensity profiles have been renormalized with their 
maximum. 



has been wide 
seismology 



observed and studied in wave physics (optics 16l. Il9l-l25l]. acoustics 26N30| 



3ll-l33l|). It manifests itself as an enhancement, by a factor 2, in the backscattered 
intensity at the vicinity of the source {i.e X = 0). Its physical origin lies in the constructive 
wave interference between reciprocal paths that have been scattered at least twice; it can 
only appear when multiple scattering occurs and the reciprocity symmetry is preserved. The 
intensity profile shown in Figj6]is a spectacular evidence of multiple scattering and shows the 
efficiency of our technique for extracting the multiple-scattering waves among a predominant 
single scattering contribution. 

Interestingly, though it is weak, the multiple scattering contribution can be taken advan- 
tage of in order to characterize the medium and determine separately the scattering losses 
and the absorption losses. When a wave propagates through a random medium, it loses 
progressively its coherence: after traveling over a distance L, only a fraction exp {—L/lext) 
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of the initial energy still propagates in coherence with the initial wave. The parameter lext, 
called the extinction mean-free path, characterizes the extinction length of the coherent part 
of the wave. It comes from two distinct phenomena (scattering and intrinsic absorption of 
the medium) which are associated to two characteristic lengths: the elastic mean-free path 
/e and the absorption mean-free path la, such that 

exp {-L/lext) = exp {-L/le) x exp {-L/Q (18) 
Experimentally, le^t can be determined b y m easurements of the ensemble-averaged field 



transmitted through a scattering layer 3J-|38|. However, this kind of experiments does not 
allow to distinguish la from le- 

We focus on the single and multiple scattering intensities obtained at the source: /"^(O, T) 
and /^^(0,T). They are plotted in FigJTl Note that the intensity of the multiple scattering 
contribution is less than 1% of the single scattering contribution. Once I^{0, T) and /^^(O, T) 
have been measured, we can fit both experimental curves with and le as independent 
adjustable parameters. To that end, we need a theoretical model describing the spatial 
and temporal evolution of the mean intensity inside the random medium. In the literature. 



the mean intensity is often assumed to obey the diffusion equation [39|. The diffusion 
approximation is simple, but only valid in the long-time limit. Since we deal with a weakly 
scattering medium, the elastic mean free path is expected to be very large compared to 
the scattering path lengths {le » cT). Thus, the diffusion approximation does not apply 
';o our problem. Instead, we used the radiative transfer equation (RTE) jo]. Paasschens 



40| proposed an exact solution of the RTE in time-domain and real space for an infinite 
2D random medium (see Appendix C), Eq.C21). Based on this theoretical work, we have 
computed the exact expression of the single and double scattering intensities, I^{0,T) and 
/(2)(0,T), as well as an approximate expression of the multiple scattering intensity I^\0, T), 
considering the medium as semi-infinite and two-dimensional. 

The choice of a 2D model is justified as follows. Experimentally, the transducers are 10 
mm in height, which is much larger than the average wavelength (0.5 mm). Moreover, a 
vertical cylindrical acoustic lens ensures that the emitted beam remains coUimated in the 
[x, z) plane. Similarly, in reception only waves propagating in the (x, z) plane are recorded by 
the transducers. Thus the single scattering problem is clearly 2D. As to multiple scattering, 
for the same reason the only paths that can generate a signal on the receiving transducers 
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FIG. 7: (a) Single-scattered intensity P{X = 0,T) versus time. Experimental measurements 
(black circles) are fitted with the theoretical curve (continuous black line) considering an extinction 
length lext = 50 mm. (b) Multiply-scattered intensity I^^ {X = 0,T) versus time. Experimental 
measurements (white squares) are compared with thoretical results for I^^{X = 0, T) (continuous 
black lines) and I^'^\X = 0,T) (dashed black lines) for different values of the mean free path, while 
keeping Ig^t = 50 mm. All intensities have been normalized by the maximum of the single scattered 
intensity I^(X = 0,T) over time. 

are those whose first and last scatterer are in the {x, z) plane. The gel sample being weakly 
scattering, I^^ is mostly dominated by double scattering (see FigjTI^b)). Thus even though 
the wave propagation in the gel sample is 3D, we have used the 2D solution for the RTE. 

The detailed calculations of J"^, I^'^^ and I^^ are shown in Appendix C. It appears that 
the single scattering intensity I^{0,T) exhibits a temporal evolution which only depends on 
the extinction length l^xt- In the case of the gel studied here, the best fit of the experimental 
results yields lext = 50 mm (FigJTl^a)). Once l^xt is known, /g and la can be determined by 
fitting l'^^{0,T) with theory (FigJT^b)) with only one adjustable parameter since 1/lext = 
l//e + l//a- The scattering gel is found to be much more absorbing than scattering : 1^ ~ 2000 
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mm, while la ~ 50 mm. FigjTl^b) also displays the theoretical evolution of the double 
scattering contribution I^'^\0,T). For /g ~ 2000 mm, J*^^^ and J*^ are nearly identical, which 
shows that the double scattering contribution clearly dominates the multiple scattering 
intensity in the gel sample. As the theoretical expression derived in Appendix C is exact for 
the measured values of and la are reliable. 
In this example, the medium was a weakly scattering gel, with a ratio I^'^ /I^ less than 
1%. Yet the separation of single and multiple scattering contributions can also be achieved 
in real scattering media for which J*'^//"^ is closer to unity, as we present in the next section. 



VI. APPLICATION TO HUMAN SOFT TISSUES 



(a) 



Breast 




FIG. 8: (a) Experimental setup used for the investigation of multiple scattering in breast tissues, 
(b) Echographic image of the breast. The grey scale is in dB. (c) Multiple scattering rate 7 as a 
function of time. 

The same kind of experiment has been performed in a biological medium for which 
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ultrasound is often used: the breast. The experimental set up is depicted in FigJHJ^a). We 
use a iV-element ultrasonic array (A^ = 125) with a 4.3 MHz central frequency and a 3.5- 
5MHz bandwidth; the array pitch 5x is 0.33 mm. The emitted signal is a 0.7-/is sinusoidal 
burst at = 4.3 MHz. The sampling frequency is 50 MHz. The experimental procedure 
is the same as in SecHIl The separation of single and multiple scattering contributions can 
be performed as in Sec JIVt but an adjustment has to be made. Indeed, in our experimental 
configuration, the array of transducers is placed in the near-field of the scattering medium 
(a = 0). Consequently, the entries of matrix K are not identically distributed: the variance 
(i.e the mean intensity /) of kij decreases significantly with the distance X = Xj — Xi between 
the source and the receiver, as shown by FigJHl This implies a different variance for each 
line of matrix A, hence modifying its theoretical distribution of singular values. The upper 
bound Amax can be computed numerically taking into account this non-uniform distribution 
of matrix elements (See Appendix D). However this is only possible at the first iteration 
{q = 1). Indeed, whereas the variance of A can be estimated, the variance distribution of 
the subspaces of A is unknown a priori. Unless we do the (strong) approximation that this 
variance is uniform, in which case EqlT3] could be applied, we cannot follow the procedure 
described in Sec JIVI Here, since A clearly has a non-uniform variance, by precaution we 
choose a different strategy: the same upper bound A^ax is considered at each iteration q 
of the single/multiple scattering separation process (see Sec JIVp . Since A^ax < Amax? this 
precaution tends to overestimate the threshold and decrease the probability of error. 

FigJH] compares the multiple scattering, single scattering and total intensity profiles, at 
a given time T. Contrary to the previous experiment in the gel sample (see FigE]), the 
spatial intensity profiles, I{X) and I^{X), are not fiat due to the near- field configuration 
of the experiment. Once again, /*^(X) exhibits a coherent backscattering peak on top of a 
fiat incoherent intensity with an enhancement factor close to 2. Interestingly, this indicates 
that K'^ is not a noise contribution, but does originate from multiple wave scattering in the 
breast tissue, even though we operate at a frequency (4.3MHz) for which human soft tissues 
are usually treated as single-scattering. 

An ultrasound image of the breast has also been obtained with the same array, using 63- 
element subapertures (FigJHJ^b)). As usual in ultrasound imaging, focusing in emission and 
reception is achieved by applying a set of 63 time delays to the signals transmitted/received 
by the array. The time delays are computed in order to focus at the desired region of 
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FIG. 9: The multiple scattering intensity I^^(X) (black dots), the single scattering intensity I^{X) 
(grey circles) and the total intensity I{X) (black squares) are plotted versus X, at time T = 35.6 
^s. The intensity profiles have been renormalized with the maximum of the total intensity. 



interest, centred at coordinates zp and Xp, assuming that the velocity of sound in soft 
tissues is known. In the case of breast, as for most soft biological tissues, it is close to 
that of water (c = 1500 m/s). Here, 2666 focal planes {zp = 8 mm to 48 mm) and 63 
values of X {xp = —316x to xp = 316x with 6x the array pitch) have been used. At 
each depth zp, a new set of time-delays is calculated; this is more demanding that classical 
ultrasound imaging techniques, which generally use the same set of time-delays as long as 
zp is within the depth of field. A line of the resulting image represents, in gray level, 
the amplitude of the total echographic signal at the focal time T = 2zp/c, once focused 
beamforming has been applied to the 63 received signals. The resulting image displays 
the reflectivity of the medium under investigation. Ultrasound images of human tissues 
usually reveal the interfaces of inner organs, and often exhibit a speckled appearance due do 
random scattering by sub-wavelength inhomogeneities (cells, fibers, tissues etc.) 4l|]. Here 
the scanned area is particularly echogene between 30-40 /is, corresponding to the depth 
range 22.5-30 mm. The typical ultrasound image of a human organ (here, the breast) is 
representative of the amplitude of backscatter at a given time which hopefully (under the 
single scattering assumption) corresponds to a given depth, but it does not allow us to 
distinguish single and multiple scattering contributions. 

However, once the matrix K has been recorded, not only can we build an echographic 
image but we can also isolate the multiple scattering contribution, and estimate the multiple 
scattering rate 7 (Eq JTO]) . 7 has been averaged over the whole frequency spectrum and is 
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displayed on FigJHJ^c) as a function of time. FigJHJ^c) complements the information brought 
by the echographic image. A relevant observation is that multiple scattering becomes pre- 
dominant from T = 46 fis, that is to say beyond a depth of 34.5 mm. It means that the 
single scattering assumption, upon which the imaging process is based, is incorrect. It does 
not mean however that the image is totally wrong; a rate 7 of 50% means that half of 
the intensity received by one individual array element comes from multiple scattering. In 
classical array imaging, each line of the picture is constructed by focused beamforming in 
emission and reception. This procedure reduces the importance of multiple scattering in the 
final image because the single scattered contributions coming from a target in the focal zone 
add up coherently whereas the contributions from multiple scattering can be expected to be 
uncorrelated. With 7 = 0.5 and assuming A^' = 63 uncorrelated array elements, the multi- 
ple scattering rate becomes 1/7/(1 — l)/N' ~ 1/63 after beamforming. This is pro bably an 
underestimation since multiple scattering signals cannot be fully uncorrelated 12|, |l5|. As 
a result, the proportion of multiple scattering in the final image around T = 50 /is is of the 
order of a few percents, which is still weak. Also note that at larger times, the technique 
we presented here would fail: after T = 50 /iS the rate of multiple scattering becomes too 
large for the SVD to extract the single scattering contribution (p = 0), at least for some 
frequencies of the spectrum. As the results are averaged over the whole frequency spectrum 
(3.5-5 MHz), the multiple scattering rates that are presented here are still meaningful, un- 
til T = 60 /is. Beyond that time, the method we presented here would be inadequate to 
separate single and multiple scattering. 

VII. CONCLUSION 

The approach we developed here can separate single- and multiple-scattered waves in 
randomly heterogeneous media. It requires an array of transmitters/receivers and takes 
advantage of the persistence of a deterministic coherence of single scattering signals along the 
antidiagonals of the inter-element matrix. Once a singular value decomposition is applied, 
the single scattering contribution (signal subspace) is separated from the multiple scattering 
contribution (noise subspace) by using a criterion based on random matrix theory. Unlike 



previous works 8 



12| this technique is particularly well-suited for weakly scattering media. 



for which single scattering dominates. In such media, the technique we presented here is 
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not intended to enhance the quahty of the ultrasound image, but rather to complement it, 
in two ways. Firstly, the experimental results indicate that this approach can be applied for 
characterization purposes: the separation of single and multiple scattering provides a way to 
measure the scattering and absorption mean-free paths independently. This idea was tested 
on a synthetic gel. Secondly, the technique was also applied in vivo to the case of breast 
imaging with ultrasonic waves around 4.3 MHz. The occurrence of multiple scattering has 
been established and its contribution to the backscattered wave-field is shown to be far from 
negligible. By measuring the relative amount of multiple scattering, the method serves as 
an experimental test for the first Born approximation (single scattering), which is usually 
made in such tissues. 
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Appendix A: 

Single and multiple scattering signals are not strictly orthogonal. Thus, S can be ex- 
pressed as the sum of the single scattering contribution plus a residual multiple scattering 
noise R: 

S = A^ + R (Al) 

The aim of this appendix is to estimate the remaining part of multiple scattering which 
corrupts the signal subspace S. In other words, we want to determine the mean intensity 
cr| of coefficients rim = (|nm|^))- 



To that aim, we apply the first order perturbation theory [42| to the autocorrelation 
matrix AA"!". Actually, AA^ can be decomposed as a sum of an unperturbed autocorrelation 
matrix A^A^^^ (linked with single scattering) and a perturbation matrix W due to multiple 
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scattering : 

AA^ = A^A^^ + W, (A2) 
with W = A^A^t + A^A^t + A^i^ (A3) 
ist order 2nd order 

As we deal with weakly scattering media, the single scattering contribution is predominant: 
as » (Jm- Thus, W can be seen as a low perturbation of the autocorrelation matrix. 
Moreover, we can neglect the second order term A'^A'^^ in Eq JA3[ W can be simplified 
into 

W ~ A^A^t + A^A^t (A4) 

Let us first focus on the reference state, i.e with no multiple scattering. In this case, we 
have 

A = A« = ^ \l^Jt^V:^^ , for aM ^ (A5) 

k=l 

where denotes the /c*^ singular value of A^. and are the singular vectors associated 
to A^. 

The perturbation theory can now be applied to AA"!". At first order, this theory states 
that only the eigenvalues (A^)^ of A A''" are perturbed by multiple scattering. The singular 
vectors remain identical to those obtained in the unperturbed case [i.e. without multiple 
scattering), hence Uk — and — V^. The p first eigenvalues (A^)^ of KA) can be 



written as 42|: 



(A,)' ^ (A^)' + U°^WU° (A6) 



(A^,)^ is the unperturbed eigenvalue (that we would obtain without multiple scattering) and 
the term U^^WU^ corresponds to the perturbation due to multiple scattering. 
Using Eqs JA6l fc IA4t one can try to express (A^)^ — (A^,)^: 

(A,)' - {Kf ^ u°^wu° 

- U^^A^A^txjo + u°^A^A^tu° 
~ A°V°^A^tuo + A°u°^A^V° 
(A,)^ - ^ 2\l^ {u°tA^V^} (AT) 
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is a random complex variable with zero mean. Let us calculate its variance: 
var [U°1'A^V°] = (u°1'A^V°V°1-A^tue' 

= U°^(A^A^t>U^ (A8) 
We assume that the matrix A^ is random, hence the ensemble average of A^A^''' yields 

(A^A^t) = Lall (A9) 
where I is the identity matrix and afj is the mean power of multiple scattering signals: 



Ml 
''Iml 



a\j. The combination of Eqs JASI fc IA9I yields: 



var 



UofAMyo 



(AlO) 



Let us calculate the variance of (A^)^ — (A^)^ (Eq JA7p : 



{[{^kf - [Kf]") ^ 4(A^)%ar 



o\2 



var 



3?iu°^A^V° 



(All) 



(A12) 



By injecting Eq lAlOl in Eq lAllt we obtain 

Using the approximation (A^)^ — (A^)^ ^ 2A^ (Afc — A^,), the last equation becomes: 

(4) - ^ (A13) 



where = A^ — A^. 

Using Eq.l5 of the manuscript and Eq JASI and the fact that Uk — UJJ and ^ V^, the 
matrix R = S — A^ can be expressed as a function of e^: 

R^X^e.U^V°^ (AM) 



k=l 



Now, we can estimate the variance ct|j of the coefficients rim- 

{\rimf) 



a 



R 



(A15) 



k=l 
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Because the singular vectors and are normalized, = 1/M and (Iffc^l^) = 1/L. 

Using Eq JA13[ we finally obtain: 

^ (A16) 



,2 



P 



^ 2M~ 

The mean power of the residual multiple scattering noise is a\ ^ V^\il (SA-f). This residual 

r| » all. 



noise is negligible compared to single scattering signals since o\» a\ 



Appendix B: 

This appendix describes the numerical process we use to determine the upper bound A^^^x 
of the singular values distribution, when short-range correlations between matrix entries 
exist. 

The first step is to estimate these correlations. This is done by measuring the correlation 
coefficients and between respectively the lines and columns of A: 

where the symbol < . > denotes an average over the variables in the subscript, i.e the 
source/receiver pairs The integer m represents the distance between sources and 

receivers, in units of the array pitch 6x. 

The second step consists in building two correlation matrices C (of size (M — g) x (M — g)) 
and D (of size Lx L) from the measured correlation coefficients F^ and F^. The coefficients 
Cij and dij of these correlation matrices are given by: 

c^J = rf_,. (B3) 
d^, = rf_,. (B4) 

Once these correlation matrices are built, the next step consists in generating a random 
matrix which displays the same short-range correlations as matrix A. 

To that aim, a random matrix P of dimension (M — g) x L is first generated numerically. 
Its coefficients are totally independent and identically distributed. Then, a matrix Q is built 
from P, such that 

Q = C^PD^ (B5) 
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The matrix Q exhibits the same correlation properties as A. The SVD of Q is then achieved 
and its first singular value Ai is renormalized according to 

Ai = , (B6) 

Y M-q l^k=l "^k 

By repeating this operation over 500 realizations, we can build a cumulative histogram of 
the first singular value Ai. The distribution function Fi of the first singular value, -Fi(A) = 
Prob{Ai < A}, can be estimated. The upper bound A^ax is then deduced from Fi. An 
acceptable probability of error 7 is first set (typically, 7 = 10~^), hence the upper bound 



iX- 

(9) 



-^max 



Ai^W = ^r'(l-7) (B7) 



Appendix C: 

In this appendix, we derive an expression for the single and multiple-scattering contribu- 
tions to the average backscattered intensity, based on radiative transfer theory. Radiative 
transfer theory describes the spatial and temporal dependence of the radiance (or specific 
intensity) P(r, t, u) in a random medium. Radiance is defined as energy fiow, propagating 
in the direction u, per unit normal area per unit solid angle dVL per unit time. It follows a 
transport (Boltzmann) equation : 

^^P(r, t, u) + u.VP(r, t, u) + Z-ip(r, t, u) = C^P(r, t) + c"^^(r, t, u), (CI) 

with S the source term and P{r,t) the intensity, defined as the angular average of the 
radiance 



Classically, the transport equation can be derived from the Bethe-Salpether equation 
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44| . neglecting all interference (coherent) effects. Here we adapt the theoretical developments 



of Paasschens 40[, which were derived for an infinite random medium, to our experimental 
configuration. The problem is solved in two dimensions. 

The typical configuration is depicted in Fig.l of the manuscript. The random medium 
{•&) is assumed to be semi-infinite, with a plane interface. It is characterized by an elastic 
mean-free path /g and an absorption length /„. The scattering events in the random medium 
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are assumed to be isotropic, hence there is no difference between the elastic mean-free path 
and the transport mean-free path /*. 

In our case, the medium is not statistically invariant under translation since it is semi- 
infinite: P depends both on the observer (r) and the source (rg). The radiative transfer 
equation is modified into: 

-|-P(r„r,t,u) + u.VP(r„r,t,u) + M(r)/-ip(r„ r, t, u) 
c at 

= M{v)l-^P{v,, r, t) + c-H{t)5{v - r,), (C3) 

where an isotropic point source is now considered: S'(r, t, u) = 5(t)(5(r — Tg). M{v) accounts 
for the semi-infinite nature size of the random medium: 

M(r') J^-"'-'^*''' (C4) 

I 0, otherwise 

This problem can be solved by considering separately the contributions to intensity from 
= 0, 1, 2, ■ ■ ■ scattering events 

oo oo 

P(rs, r, t, u) = ^ Pjv(r„ r, t, u), P(r„ r, t) = Pn{v,, r, t) (C5) 

N=0 N=0 

The partial intensities Pat satisfy 

(^A + u.V + M(r)e,) Po(r„ r, t, u) = c-i5(t)5(r) (C6) 

(^-^ + u.V + M(r)C,i^P^(r„r,t,u) = M{rX'Pj,^,{r,,r,t), for N > (C7) 

Let us first investigate the ballistic intensity Pq. Following Paaschens 40|, the differential 
operators on the left hand side of Eq JC6l can be integrated, to yield 

1 /"°° 

Po(rs,r,t,u) = - / c/roe-^("-"^+"°")/'-*5(r - rou - rg)5(t - ro/c) 
c Jo 

^ g-/(r.,r.+ctu)/«e.t^(j. _ _ j.^) 

-f(rs,r)/lext 

= ^5(|r - r,| - ct)6 (u - u,,r) (C8) 

r\r-rs\ 

where /(rs,r) = / dr' M {r^ + r'u^^^) (C9) 
Jo 
r - 



and u 



s,r 



|r - 

The Dirac distribution 5(|r — rs| — ct) describes the propagation of the energy pulse at the 
finite speed c. |r — rs|~^ accounts for the geometrical spreading of the ballistic radiance in 
two dimensions, /(fg, r) is the ballistic path length within the scattering medium: 
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• If Ts, r e {-&), then /(rg, r) = |r - rs|. 

• Ifrs,r^ (^?),then/(rs,r) = 0. 

• In other configurations (see Fig JTOl) . the balhstic radiance propagates partly in free 
space, partly in the random medium: < /(rg, r) < |r — rs|. 

The angular average of Eq lUSI yields the balhstic intensity Po{rs,r,t): 

p-f{rs,r)/lext 

Po(rs, r, t) = — -6{\v - r,| - ct) (CIO) 



27r r — r 



Now let us consider the scattered contributions (A^ > 1). Similarly, the differential 
operators on the left hand side of Eq JC7l can be integrated, to yield 



-P/v(rs, I"' 



t, u) = / drol-^M{v^ + rou)e-^('"=''"=+'"o")/''=-'Pw_i(rs + TqU, r, t - ro/c) 



P^(r„r,t) = - / / dh^M{T^)— ^p^_,(ri,r,t-|ri-r,|/c) (Cll) 

k J J 27r|ri-rs| 

In this expression, the first scattering takes place in ri, somewhere inside {'&). 

g-/(rs,ri)//e3:t^^27r|ri — Fsl) corrcspouds to the amplitude of the ballistic intensity at ri. The 

ratio d?vi/l(. can be interpreted as a scattering cross section (actually a length, in 2D) of the 

infinitesimal surface element d?v-i. The intensity propagation from ri to r and the associated 

— 1 scattering events are taken into account by Pn~i- Hence the partial intensities can 
be obtained recursively. 

From the expression of Pq (Eq JC10(l . we infer the single-scattering term Pi: 



I r r g-/(rs,ri)/«ea:t 

Pi(rs,r,t) = -r rf2riM(ri)— -Po(ri, r, t - In - r^l/c) 

k J J 27r|ri-rs| 

leJ J 27r|ri-r,| 27r|r-ri| 

This expression can be used to calculate the single scattering contribution for comparison 
with the experimental measurements. The transducers have a finite size fe, hence a directivity 
pattern ViO): 

V{e) = ^i^W^tang/c) ^^^3^ 
Trbfctan6/c 

where fc is the central frequency. The same transducer is used as a source and a receiver. 
For convenience, it is placed at the origin: rg = r = 0. Hence from Eq lC12l we have 

I\X = 0,t) = P,{O,O,t) = ^^ 11 rfri|P(^i)l -2 Si2n-ct) (C14) 
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Using polar coordinates ri = (ri, 9i) and with /(O, ri) = ri — -^^^ (see Fig JTOj) . we obtain: 

with 9max = ^rccos (^) and 6 the Heaviside function. Eq JClSI is the simplest analytical 
form that we can obtain for the single- scattering intensity in our configuration. I^{t) can 
only be computed numerically (see Fig. 6 (a) of the manuscript). Once it is normalized by 
its maximum, the evolution of the single scattering intensity with time only depends on the 
extinction length l^xt- 




FIG. 10: Typical configuration used for the analytical calculation of the single scattering intensity 
I^{X = 0,t) 

We now turn to higher orders of scattering, starting with N = 2. Combining Eq JClll and 
IC12I , we obtain: 

1 r r p-f{rs,ri)/lext r r 

P2(rs,r,t) = j dh^M{r^) ^^^^^ _ ^ ^ j j dh^Miv^) 

g-/(r2,r)/«erfg-/(ri,r2)/«eit 

X -7r\ k^n p5(|ri - Fsl + |r2 - ril + |r- ral - ct) (C16) 

2tt\y2 — ri|27r|r — 

Rigorously, we would have to calculate all partial intensities of order N to obtain a theoretical 
expression of the multiple scattering intensity P^^(rs,r,t) = ^^^2 -Pi(rs5 r, t). This would 
be particularly tedious, time-consuming numerically, and useless in our experimental case. 
A simplified way is proposed to approach the multiple scattering intensity. Summation of 



28 



Eqs lCTOllCTTl over all N results in: 



2 M(ri)e--'^(''^''"i)/^-' 



P(rs,r,t) = Po(rs,r,t)+ / / d'r^ '^^^ P(ri, r, t - In - r^l/c) (C17) 



2 M(ri)e--^('"-''i)/'-* 



Po(rs,r,t)+ / / rf^ri ^ Po(ri, r, t - In - r,|/c) 

27r/e ri - Fs 



A(rs,r,i) 



2 M(n)e-^('"=''"^)/'=^* 



rf'ri^-^^^ [P(n, r, t - In - rsl/c) - Po(n, r, t - |n - rs|/c|pl8) 



The expression of P*^ (Eq JC18() can be rewritten using the spatial reciprocity theorem. The 
vector positions r and ri can be exchanged in P(ri,r,t — |ri — rs|/c). P^^ then becomes 

P^^(n,r,t) = - / / dh, Mir,"" 



k J J 27r|ri-rs| 

X [P(r, n, t - |ri - r,|/c) - Po(r, n, t - - n|/c)](C19) 

Replacing P(r, ri, t — |ri — rs|/c) by its expression in Eq JClTI finally leads to: 
P^(n,r,t) = ^ / / rfVM(n) 



IIJ J 27r|ri - n 

Q--f{r,r2)/lcxt 



X / / d'r^M{T^)— ^P(ri,r2,t-(|n-rs| + |r-r2|)/c) (C20) 

In this expression of P^^, ri and T2 can be interpreted as the points where the first and 
the last scattering events take place. The propagation of the intensity between ri and T2 is 
described by P. Eq JC20l is still implicit since P^'^ is actually contained in P. An approximate 
solution for E( jC20l can be obtained by taking for P the solution of the RTE for the case of 



an infinite random medium 



40| 



Poo(ri,r2,t) = Pocir = \r2- ri\,t) 

5{ct - r) 



27rr 



ballistic term 

p-ct/lext / j,2 \ -1/2 / ^/„2f2 _ 

^^V-^) exp ^^^_|e(c«-.) (C21, 



scattering term 

Replacing P by Pqo in Eq JC20l allows to approach the multiple scattering intensity P^^: 

1 r r p-f{rs,ri)/lext 

P^^(n,r,t) ^ j^JJ rf^riM(n) 2^1^^ _^ I 

r r p-f{r,r-2)/lext 

X / / d''v2M{v2)— -Poo(ri,r2,t-(|n-rs| + |r-r2|)/c)(C22) 

J J 27r|r-r2| 
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One can object that it would have been simpler to replace P by Poo in Eq JClSI Yet this 
would amount to considering scattering paths whose last scatterer is outside of which is 
unphysical. As a consequence, the double scattering contribution would have been largely 
overestimated. On the contrary, Eq JC22l only takes into account scattering paths for which 
the first and last scatterers are within the random medium {'&). Therefore Eq JC22l exactly 
accounts for the double scattering contribution. Yet it overestimates higher orders of scat- 
tering. 

The directivity patterns of the transmitting/receiving elements as well as the distance 
a between the array and the sample can be taken into account in the calculation of the 
multiple scattering intensity, in a similar manner as what was done for the single scattering 
contribution. At exact backscattering {X = 0), because of the coherent backscattering effect 



231], the multiple scattered intensity is twice that predicted from RTE. Hence we have, for 



the double scattering intensity : 
= 0,t) = 2P2(0,0,t) 



27rri 



dW 



\V( 



-(/(r2,0)+/(ri,r2))/«erf 



'{^) 47r^|r2 - ri|r2 

with /(O, Ti) = Ti — a/ cos 9i and /(ri, = |r2 — ri|. 



-5(ri + |r2-ri|+r2-ctOC23) 




FIG. 11: Typical configuration used for the analytical calculation of the double scattering intensity 

J(2)(0,t) 

The Dirac function which appears upon the integrand of Eq JC23l implies that the integral 
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over ri has to be performed over the intersection between the disk V[0, ct/2] (whose center 
is and radius is ct/2) and the random medium (^9) (the grey cross hatched surface in 
Fig JTT]) . The integration over r2 is performed along the intersection between the elhpse 
£^[0,ri, {ct — ri)/2] (whose foci are and ri, and semi major axis is (ct — Ti)/2) and the 
random medium {^) (see Fig JTT]) . Using polar coordinates, Eq JC23l becomes 

-ct/lext r r 

/(2)(X = 0,t) = / / drirf^i |P(ei)|'e"/(^-*™^^^) 



Eq JC24l is the simplest analytical form for the double-scattering intensity in our configura- 
tion. J*^^) can be computed numerically (see Fig. 6(b) of the manuscript). 

Finally, we consider the approximate solution for the multiple scattering intensity 
(Eq JC22( ). and take into account the coherent backscattering phenomenon, as well as the 
directivity of the sources/receivers to obtain 

IIJ J{^{~ 27rri 



Using polar coordinates, Eq lC25l becomes 



, p-/{0,r2)//e.i 

X // dh2\V{92)\' — Poo(|r2-ri|,t-(ri + r2)/c) (C25) 



2tt L .1 Jvf]{&) 



X 



/ / dr2de2 \V{e2)f e-("^-'^/™^'^)/^^-Poo(|r2 - ri|, t - (n + r2)/cj(C26) 

Unlike Eqs JClSl and \C2A\ Eq JC26l is an approximate solution of the radiative transfer equa- 
tion. It is the simplest analytical expression that we can obtain. It allows to compute the 
time-evolution of /^^ numerically (see Fig.6(b)). Note that J*^(t) does not only depend on 
lext, but also on the relative weight of the diffusion (/g) and the absorption (/„) losses. This 
dependence on 1^ and is contained implicitly in the intensity propagator Poo(see Eq jC2ip . 



Appendix D: 

This appendix describes the numerical process we use to determine the upper bound 
Amax of the singular values distribution, when the variance of the element aij is not the 
same for every line of the matrix A. 
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The first step is to estimate this variance. This is done by measuring the mean intensity 
li of Qij along each hne z of A: 

h = {\a,f)^ (Dl) 

where the symbol < . > denotes an average over j. 

The second step consists in building a diagonal matrix S of size M x M, whose elements 
are: 

<Tii = \fTi (D2) 

The next step is the generation of a random matrix Q with the same variance as A. 
To that aim, a random matrix P of dimension M x L is first generated numerically. Its 
coefficients are independent and identically distributed. Then Q is built from P according 
to : 

Q = SC^PD^ (D3) 

This ensures that the g^j have the same variance and correlation properties than the 
aij. Once the SVD of Q is performed, its first singular value Ai is renormalized according 
to Eq.lO. The same operation is repeated for 500 realisations, which yields a histogram 
of Ai and thereby an estimate of the distribution function F\ of the first singular value, 
Fi(A) = Prob{Ai < A}, can be estimated. The upper bound A^ax is then deduced from F\ 
(see Appendix [Bl) . 
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